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Abstract 



It is well-known that spatial averaging can be realized (in space or 
Y \ frequency domain) using algorithms whose complexity does not scale 

with the size or shape of the filter. These fast algorithms are generally 
referred to as constant- time or O(l) algorithms in the image processing 
literature. Along with the spatial filter, the edge-preserving bilateral 
filter [12] involves an additional range kernel. This is used to restrict 
the averaging to those neighborhood pixels whose intensity are similar 
or close to that of the pixel of interest. The range kernel operates by 
acting on the pixel intensities. This makes the averaging process non- 
CNJ linear and computationally intensive, especially when the spatial filter 

is large. In this paper, we show how the O(l) averaging algorithms 
can be leveraged for realizing the bilateral filter in constant-time, by 
using trigonometric range kernels. This is done by generalizing the 
idea in |10| of using polynomial kernels. The class of trigonometric 
I kernels turns out to be sufficiently rich, allowing for the approximation 

of the standard Gaussian bilateral filter. The attractive feature of our 
k> approach is that, for a fixed number of terms, the quality of approxi- 

j_j mation achieved using trigonometric kernels is much superior to that 

obtained in [10 1 using polynomials. 
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1 Introduction 



The bilateral filtering of an image f(x) in the general setting is given by 

f(x) = r]- 1 J w(x, y) <p(f(x), f(y)) f(y) dy 

where 

77 = J w(x, y) <p(f(x), f{y)) dy. 

In this formula, w(x,y) measures the geometric proximity between the pixel 
of interest x and a nearby pixel y. Its role is to localize the averaging to 
a neighborhood of x. On the other hand, the function 4>{u, v) measures 
the similarity between the intensity of the pixel of interest f(x) and its 
neighbor f(y). The normalizing factor n is used to preserve constants, and 
in particular the local mean. 

In this paper, we consider the so-called unbiased form of the bilateral 
filter [12J, where w(x, y) is translation-invariant, that is, w(x, y) = w(x — y), 
and where the range filter is symmetric and depends on the difference of 
intensity, 4>(f(x), f(y)) = 4>(f(x) — f{y)). In this case, the filter is given by 

f(x) = tT 1 / w(y)cP(f(x -y)- f(x))f(x - y) dy (1) 
Jn 

where 

n = w{y)4>{f(x -y)- fix)) dy. (2) 
Jn 

We call w(x) the spatial kernel, and (f>(s) the range kernel. The local support f2 
of the spatial kernel specifies the neighborhood over which the averaging 
takes place. A popular form of the bilateral filter is one where both w{x) 
and ct>(s) are Gaussian |T2HIDH2J|TS|. 

The edge-preserving bilateral filter was originally introduced by Tomasi 
et al. in fT2| as a simple, non-iterative alternative to anisotropic diffusion 
18 J. This was motivated by the observation that while standard spatial 
averaging performs well in regions with homogenous intensities, it tends 
to performs poorly in the vicinity of sharp transitions, such as edges. For 
the bilateral filter in ([!), the difference f(x — y) — f{x) is close to zero in 
homogenous regions, and hence 4>{f{x — y) — f{x)) rj 1. In this case, ([TJ 
simply results in the averaging of pixels in the neighborhood of the pixel of 
interest. On the other hand, if the pixel of interest a; is in the vicinity of an 
edge, 4>(f(x — y) — f{x)) is large when x — y belongs to the same side of 
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the edge as x, and is small when x — y is on the other side of the edge. As 
a result, the averaging is restricted to neighborhood pixels that are on the 
same side of the edge as the pixel of interest. This is the basic idea which 
allows one to perform smoothing while preserving edges at the same time. 
Since its inception, the bilateral filter has found widespread use in several 
image processing, computer graphics, and computer vision applications. 
This includes denoising JT], video abstraction [14], demosaicing JIT], optical- 
flow estimation JT5], and stereo matching [17], to name a few. More recently, 
the bilateral filter has been extended by Baudes et al. |2] to realize the 
popular non-local neighborhood filter, where the similarity between pixels 
is measured using patches centered around the pixels. 

The direct implementation of ([lj turns out to be rather computationally 
intensive for real time applications. Several efficient numerical schemes 
have been proposed in the past for implementing the filter in real time, 
even at video rates [5l [13j |9j [7|- These algorithms (with the exception of 
13), however, do not scale well with the size of the spatial kernel, and this 
limits their usage in high resolution applications. A significant advance 
was obtained when Porikli fj"0| proposed a constant-time implementation 
of the bilateral filter (for arbitrary spatial kernels) using polynomial range 
kernels. The O(l) algorithm was also extended to include Gaussian <p(s) by 
locally approximating it using polynomials. More recently, Yang et al. [fl"6| 
have proposed a 0(1) algorithm for arbitrary range and spatial kernels by 
extending the bilateral filtering method of Durand et al. [5J. Their algorithm 
is based on a piecewise-linear approximation of the bilateral filter obtained 
by quantizing <fi(s). 

In this paper, we extend the 0(1) algorithm of Porikli to provide an exact 
implementation of the bilateral filter, using trigonometric range kernels. Our 
main observation that trigonometric functions share a common property 
of polynomials which allows one to "linearize" the otherwise non-linear 
bilateral filter. The common property is that the translate of a polynomial 
(resp. trigonometric function) is again a polynomial (resp. trigonometric 
function), and importantly, of the same degree. By fixing </>(s) to be a 
trigonometric function, we show how this self-shiftable property can be 
used to (locally) linearize the bilateral filter. This is the crux of the idea that 
was used for deriving the 0(1) algorithm for polynomial 4>(s) in [10|. 
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2 Constant-time bilateral filter 



2.1 The main idea 

It is the presence of the term 4>(f(x — y) — f(x)) in ([I) that makes the filter 
non-linear. In the absence of this term, that is, when <fi(s) is constant, the 
filter is simply given by the averaging 

J(x) = / w(y)f(x - y) dy, (3) 
Jn 

where we assume w(x) to have a total mass of unity. It is well-known 
that ^ can be implemented in constant-time, irrespective of the size and 
shape of the filter, using the convolution-multiplication property of the 
(fast) Fourier transform. The number of computations required per pixel, 
however, depends on the size of the image in this case H51 . On the other 
hand, it is known that <|3j can be realized at the cost of a constant number 
of operations per pixel (independent of the size of the image and the filter) 
using recursive algorithms. These 0(1) recursive algorithms are based on 
specialized kernels, such as the box and the hat function |6ll4llT9|, and the 
more general class of Gaussian-like box splines 1151 . 

Our present idea is to leverage these fast averaging algorithms by ex- 
pressing ([IJ in terms of j3|, where the averaging is performed on the image 
and its simple pointwise transforms. Our observation is that we can do so if 
the range kernel is of the form 

4>{s) = cos(7s) (-T <s<T). (4) 

By plugging Q into Q, we can write the integral as 

cos(7/(s)) / w(y)cos(jf(x-y))f(x-y)dy+sm(jf(x)) / w(y) sm(-jf(x-y))f(x-y) dy. 
Jn Jn 

This is clearly seen to be the linear combination of two spatial averages, 
performed on the images cos^ f (x)) f (x) and sin(7/(a;))/(a5). Similarly, we 
can write the integral in §2§ as 



cos(7/(£c)) / w(y) cos(j f(x-y)) dy+sm(>yf(x)) / w(y) sin(j f(x-y)) dy. 
Jn Jn 

In this case, the averaging is on the images cos(7/(a;)) and sin(7/(a;)). This 
is the trick that allows us to express Q in terms of linear convolution filters 
applied to pointwise transforms of the image. 
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Note that the domain of 4>(s) is [— T, T] in We assume here (without 
loss of generality) that the dynamic range of the image is within [0, T]. The 
maximum of \f(x) — f(y) | over all x and y such that x — y £ fHs within T in 
this case. Therefore, by letting 7 = tt/2T, we can guarantee the argument 7s 
of the cosine function to be within the range [— tt/2, 7r/2]. The crucial point 
here is that the cosine function is oscillating and can assume negative values 
over (—00, 00). However, its restriction over the half-period [—tt/2, tt/2] has 
two essential properties of a range kernel — it is non-negative and has a 
bump shape (cf. the outermost curve in Figure [TJ. Note that, in practice, the 
bound on the local variations of intensity could be much lower than T. 

2.2 General trigonometric kernels 

The above idea can easily be extended to more general trigonometric func- 
tions of the form cf>(s) = ao + a\ cos(7s) + • • • + cos(N^fs). This is most 
conveniently done by writing <p(s) in terms of complex exponentials, namely 
as 

4>{ s ) = ^2 Cn exp n ^ s ) ■ ( 5 ) 

\n\<N 

The coefficients c n must be real and symmetric, since (p(s) is real and sym- 
metric. Now, using the addition-multiplication property of exponentials, 
we can write 

c/)(f(x - y) - f{x)) = d n{x) exp (jnjf(x - y)) 

\n\<N 

where d n (x) = c n exp ( — jnjf(x)). Plugging this into Q, we immediately 
see that 

7, N J2\n\<N d n {x) g n {x) 

j(x) = = (6) 

zZ\n\<N d n (x) h n (x) 

where h n (x) = exp (jnjf(x)), and g n (x) = f(x)h n (x). We refer to h n (x) 
and g n (x) as the auxiliary images, and N as the degree of the kernel. 

The above analysis gives us the following 0(1) algorithm for the bilateral 
filter: We first set up the auxiliary images and the coefficients d n [x) from 
the input image. We then average each of the auxiliary images using a 0(1) 
algorithm (this can be done in parallel). The samples of the filtered image 
is then given by the simple sum and division in j6|. In particular, for an 
image of size M x M, we can compute the spatial averages for any arbitrary 
w(x) at the cost of 0(M 2 log 2 M) operations using the Fourier transform. As 
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Figure 1: The family of raised cosines g(s) = [cos(7s)] over the dynamic 
range — T < s < T as N goes from 1 to 5 (outer to inner curves). We set 
T = 255 corresponding to the maximum dynamic range of a grayscale 
image, and 7 = ir/ 2T. They satisfy the two essential properties required 
to qualify as a valid range kernel of the bilateral filter — non-negativity and 
monotonicity (decay). Moreover, they have the remarkable property that 
they converge to a Gaussian (after appropriate normalization) as N gets 
large; see 0. 
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Figure 2: Approximation of the Gaussian exp(— x 2 /2a 2 ) (dashed black curve) 
over the interval [—255, 255] using the Taylor polynomial (solid red curve) 
and the raised cosine (solid blue curve). We set a = 80, and use N = 4 for the 
raised cosine in (|7|). The raised cosine is of the form ao+ai cos(2#)+a2 cos(4#) 
in this case. We use a 3-term Taylor polynomial of the form bo + b\x 2 + 62X . 
It is clear that the raised cosine offers a much better approximation than its 
polynomial counterpart. In particular, note how the polynomial blows up 
beyond \x\> 100. 

mentioned earlier, this can further be reduced to a total of 0(M 2 ) operations 
using specialized spatial kernels (HIHEI. 

2.3 Raised cosines 

We now address the fact that <p(s) must have some additional properties to 
qualify as a valid range kernel (besides being symmetric). Namely, 4>(s) must 
be non-negative, and must be monotonic in that <fi(si) < 4>{s2) whenever 
|si| > |s2 1 . In particular, it must have a peak at the origin. This ensures that 
large differences in intensity gets more penalized than small differences, and 
that ([TJ behaves purely as a spatial filter in a region having uniform intensity. 
Moreover, one must also have some control on the variance (effective width) 
of 4>(s). We now address these design problems in order. 
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The properties of symmetry, non-negativity and monotonicity are simul- 
taneously enjoyed by the family of raised cosines of the form 

<j>(s) = [cos(-fs)] N (-T<s<T). 

Writing cos 6 = [e? e + e~^ e )/2, and applying the binomial theorem, we see 
that 

N . . 

4>(s) = £ 2- N r) exp (i(2n - N) 7 s). 

n=0 \ n ' 

This expresses the raised cosines as in ||5j, though we have used a slightly 
different summation. Since 4>{s) has a total of (N + 1) terms, this gives a total 
of 2{N + 1) auxiliary images in d6). The central term n = N/2 is constant 
when N is even, and we have one less auxiliary image to process in this 
case. 



2.4 Approximation of Gaussian kernels 

Figure [T] shows the raised cosines of degree N = 1 to N = 5. It is seen 
that <p(s) become more Gaussian-like over the half-period [—it, it] with the 
increase in N. The fact, however, is that <p(s) converges pointwise to zero at 
all points as N gets large, excepting for the node points 0, ±7r, ±2tt, . . .. This 
problem can nevertheless be addressed by suitably scaling the raised coinse. 
The precise result is given by the following pointwise convergence: 



lim 

N — >oo 
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Proof. Note that Taylor's theorem with remainder tells us that if f(x) is 
sufficiently smooth, then f(x) = JX^ 1 x fe / (fc) (0) /k\ + x m f^ {6) /m\, where 
9 is some number between and x. Applied to the cosine function, we have 
cos(x) = 1 — x 2 /2 + x 4 cos 0/24. In other words, cos(x) = 1 — x 2 /2 + r(x), 
where \r(x)\ < x 4 (we write f(x) < g(x) to signify that f(x) < Cg{x) for 
some absolute constant C, where C is independent of x). Using this estimate, 
along with the binomial theorem, we can write 
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where |r(s, N)\ < s A /N 2 . We are almost done since it is well-known that 
(1 + x/N) N approaches exp(x) as N gets large. To establish 0, all we need 
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to show is that, for any fixed s, the residual terms can be made negiligibly 
small simply by setting N large. 

Now note that if \s\ < N 1 / 2 , then the magnitude of (1 - ^ 2 s 2 /2N) is 
within unity and, on the other hand, s 4 /N < 1 when |s| < iV 1 / 4 . Thus, 
given any fixed s, we set N to be large enough so that s satisfies the above 
bounds. Then, following the trivial inequality (*P\ < N k , we see that the 
modulus of the residual is 

N / 9 4 \ k /^\ N 1 
< vW — ] <N — ) <- 
\ N *J ~ \ N J ~ Ar ' 

provided that |s| < Ln = (N 1 - 2 ^) 1 ^. This can clearly be achieved by 
increasing N, since Ljv is monotonic in N. □ 

We have seen that raised cosines of sufficiently large order provide ar- 
bitrarily close approximations of the Gaussian. The crucial feature about 
^ is that the rate of convergence is much faster than that of Taylor polyno- 
mials, which were used to approximate the Gaussian range kernel in IfTUI . 
In particular, we can obtain an approximation comparable to that achieved 
using polynomials using fewer number of terms. This is important from 
the practical standpoint. In Figure |2j we consider the target Gaussian kernel 
exp(— s 2 /2a 2 ), where a = 80. We approximate this using the raised cosine 
of degree 4, which has 3 terms. We also plot the polynomial corresponding 
to the 3-term Taylor expansion of the Gaussian, which is used in for approx- 
imating the Gaussian in [10] . It is clear that the approximation quality of 
the raised cosine is superior to that offered by a Taylor polynomial having 
equal number of terms. In particular, note that the Taylor approximation 
does not automatically offer the crucial monotonic property. 

Table 1: iVo is the minimum degree of the raised cosine required to approx- 
imate a Gaussian of standard deviation a on the interval [—255, 255]. The 
estimate [(7c) -2 ] is also shown. 
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2.5 Control of the width of range kernel 

The approximation in Q also suggests a means of controlling the variance of 
the raised cosine, namely, by controlling the variance of the target Gaussian. 
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The target Gaussian (with normalization) has a fixed variance of j~ 2 . This 
can be increased simply by rescaling the argument of the cosine in §7§ by 
some p > 1. In particular, for sufficiently large N, 



cos 



7 s 



pVN 



' s 2 



eX P ( -^2^2 J • ( 8 ) 



The variance of the target Gaussian (again with normalization) has now 
increased to p 2 j~ 2 - A fairly accurate estimate of the variance of the raised 
cosine is therefore a 2 m p 2 j~ 2 - In particular, we can increase the variance 
simply by setting p = for all a > 7 -2 , provided iV is large enough. 

Bringing down the variance below j~ 2 , on the other hand, is more subtle. 
This cannot be achieved simply by rescaling with p < 1 on account of the 
oscillatory nature of the cosine. For instance, setting p < 1 can cause (f>(s) to 
become non-negative, or loose its monotonicity. The only way of doing so is 
by increasing the degree of the cosine (cf. Figure [lj. In particular, N must 
be large enough so that the argument of cos(-) is within [— ir/2, it/ 2] for all 
-T < s<T. This is the case if 



iV > P ~ (7cr) . 

In other words, to approximate a Gaussian having a small variance a, N 
must at least be as large as Nq ~ (70-) -2 . The bound is quite tight for large 
a, but is loose when a is small. We empirically determined No for certain 
values of a for the case T = 255, some of which are given in Table [l] It 
turned out to be much lower than the estimate (ja)^ 2 when a is small. For 
a fixed setting of T (e.g., for grayscale images), this suggests the use of a 
lookup table for determining iVo for small a on-the-fly 

The above analysis leads us to an 0(1) algorithm for approximating 
the Gaussian bilateral filtering, where both the spatial and range filters are 
Gaussians. The steps are summarized in Algorithm [l] 



3 Experiments 

We implemented the proposed algorithm for Gaussian bilateral filtering in 
Java on a Mac OS X 2x Quad core 2.66 GHz machine, as an ImageJ plugin. 
We used multi-threading for computing the spatial averages of the auxiliary 
images in parallel. A recursive 0{\) algorithm was used for implementing 
the Gaussian filter in space domain [18|. The average times required for 
processing a 720 x 540 grayscale image using our algorithm are shown in 
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Algorithm 1 Fast 0(1) bilateral filtering for the Gaussian kernel 

Input: Image f(x), dynamic range [— T, T], a 2 and a 2 , for the spatial and 
range filters. 

1. Set 7 = tt/2T, and p = 7<7 r . 

2. If a r > 7~ 2 , pick any large N. Else, set N = (ja r )~ 2 , or use a look-up 
table to fix N. 

3. Set h n (x) = exp (jj{2n - N)f(x)/py/Nj and g n (x) = f(x)h n (x), and 

the coefficients d n {x) = 2~ N (*) exp (-n(2n - N) f \x) / py/N) . 

4. Use an 0(1) algorithm to filter h n {x) and g n (x) with a Gaussian of 
variance a 2 to get h n (x) and g~^{x). 

5. Set f(x) as the ratio of J2 n =o d n (x)g^(x) and Y^n=o d n (x)h n (x). 
Return: Filtered image f(x). 



Table |2] We repeated the experiment for different variances of the Gaussian 
range kernel, and at different spatial variances. As seen from the table, 
the processing time is quite fast compared to a direct implementation of 
the bilateral filter, which requires considerable time depending on the size 
of the spatial filter. For instance, a direct implementation of the filter on 
a 512 x 512 image required 4 seconds for a s as low as 3 on our machine 
(using discretized Gaussians supported on [—3a, 3a] 2 ), and this climbed up 
to almost 10 seconds for a s = 10. As is seen from Table|2j the processing time 
of our algorithm, however, suddenly shoots up for narrow Gaussians with 
o r < 15. This is due to the large N required to approximate the Gaussian in 
this regime (cf. Table [ll. We have figured out an approximation scheme for 
further accelerating the processing for very small oy, without appreciably 
degarding the final output. Discussion of this method is however beyond 
the present scope of the paper. 

We next tried a visual comparisonof the ouput of our algorithm with the 
algorithm in [10 j. In Figure |3j we compare the outputs of the two algorithms 
with the direct implementation, on a natural grayscale image. As is clearly 
seen from the processed images, our result resembles the exact output very 
closely. The result obtained using the polynomial kernel, on the other 
hand, shows strange artifacts. The difference is also clear from the standard 
deviation of the error between the exact output and the approximations. We 
note, however, that the execution time of the polynomial method is slightly 
lower than that of our method, since it requires half the number of auxiliary 
images for a given degree. 

We also tested our implementation of the Gaussian bilateral filter on 
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color (RGB) images. We tried a naive processing, where each of the three 
color channels were processed independently. The results on a couple of 
images are shown in Figure |4] The Java source code can be downloaded 
from the web at http:/ /bigwww.epfl.ch/algorithms/bilateral-filter 



Table 2: The time in milliseconds required for processing a grayscale image 
of size 720 x 540 pixels using our algorithm. The processing was done on a 
Mac OS X, 2x quad core 2.66 GHz machine, using multithreading. 
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4 Discussion 

We presented a general method of computing the bilateral filter in constant- 
time using trigonometric range kernels. Within this framework, we showed 
how feasible range kernels can be realized using the family of raised cosines. 
The highlights of our approach are the following: 

• Accuracy. Our method is exact, at least for the family of raised cosines. It 
does not require the quantization of the range kernel, as is the case in Il5lll~6l. 
Moreover, note that the auxiliary images in Q have the same dynamic range 
as the input image irrespective of the degree N. This is unlike the situation 
in [10], where the dynamic range of the auxiliary images grow exponentially 
with the N. This makes the computations susceptible to numerical errors 
for large N. 

• Speed. Besides having 0(1) complexity, our algorithm can also be imple- 
mented in parallel. This allows us to further accelerate its speed. 

• Approximation property. Trigonometric functions yield better (local) 
approximation of Gaussians than polynomials. In particular, we showed 
that by using a particular class of raised cosines, we can obtain much better 
approximations of the Gaussian range kernel than that offered by the Taylor 
polynomials in [To). The final output is artifact-free and resembles the true 
output very closely. The only flip side of our approach (this is also the case 
with [10], as noted in [16J) is that a large number of terms are required to 
approximate very narrow Gaussians over large intervals. 

• Space-variant extension. The spatial kernel in Q can be changed from 
point-to-point within the image to control the amount of smoothing (partic- 
ularly in homogenous regions), while the range kernel is kept fixed. Thanks 
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(a) (b) 




(c) (d) 

Figure 3: Comparison of various implementations of the Gaussian bilateral 
filter on the grayscale image Isha of size 600 x 512. The filter settings are 
a s = 15 and a r = 80. (a) Original image; (b) Direct implementation of 
the bilateral filter; (c) Output obtained using polynomial kernel [10J; and 
(d) Output of our algorithm. Note the strange artifacts in (c), particularly 
around the right eye (see zoomed insets). This is on account of the distortion 
caused by the polynomial approximation shown in Figure [2] The standard 
deviation of the error between (b) and (c) is 6.5, while that between (b) and 
(d) is 1.2. 
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Figure 4: Results on the color images Greekdome and Tulip, using our im- 
plementation of the Gaussian bilateral filter. The original image is on the 
left, and the processed image is on the right. In either case, the red, green, 
and blue channels were processed independently. We used a s = 10 and 
o r = 20 for Greekdome, and a s = 20 and a r = 60 for Tulip. (Images courtsey 
of Sylvain Paris and Fredo Durand). 
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to this can be done simply by computing the space-variant averages of 
each auxiliary image. The good news is that this can also be realized for 
a M x M image at the cost of 0(M 2 ) operations, using particular spatial 
kernels. This includes the two-dimensional box and hat filter [6. 4j, and the 
more general class of Gaussian-like box splines in [3j. 
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